---
title: "SEM: Rights Awareness, Bullying Experience, and Mental Well-Being"
output: pdf_document
---

```{r setup, include=FALSE}
library(tidyverse)
library(haven)
library(lavaan)
library(semPlot)
library(kableExtra)
library(psych)
library(semTools)
```

## 1. Load Data

```{r load-data}
klt <- read_dta("klt2024.dta")
```

## 2. Specify SEM Model

```{r model-specification}
model <- '
  # Latent variables
  RightsAwareness =~ CREDREL + CREDMEC + CREDAGE + CREDGEND + CREDDIS + CREDPOL
  BullyingExperience =~ BULLFMEV1 + BULLFMEV2 + BULLFMEV5 + BULLFMEV6
  WellBeing =~ KY10IN01 + KY10IN02 + KY10IN03 + KY10IN04 + KY10IN05 +
               KY10IN06 + KY10IN07 + KY10IN08 + KY10IN09 + KY10IN10

  # Structural paths
  WellBeing ~ RightsAwareness + BullyingExperience
'
```

## 3. Fit SEM Model

```{r model-fit}
fit <- sem(model, data = klt, missing = "fiml")
summary(fit, fit.measures = TRUE, standardized = TRUE)
```

## 4. Enhanced Path Diagram (Embedded + Exported)

```{r path-diagram, fig.align='center', fig.height=8, fig.width=14}
# Export PNG version
png("path_diagram.png", width = 2400, height = 1400, res = 220)
semPaths(fit,
         what = "std",
         layout = "tree",
         edge.label.cex = 0.9,
         label.font = 1,
         sizeMan = 6,
         sizeLat = 8,
         residuals = FALSE,
         intercepts = FALSE,
         thresholds = FALSE,
         nCharNodes = 0,
         mar = c(6, 6, 6, 6),
         color = list(lat = "lightblue", man = "white"),
         edge.color = "black",
         edge.width = 1.5,
         title = FALSE)
dev.off()

# Embed same diagram in PDF output
semPaths(fit,
         what = "std",
         layout = "tree",
         edge.label.cex = 0.9,
         label.font = 1,
         sizeMan = 6,
         sizeLat = 8,
         residuals = FALSE,
         intercepts = FALSE,
         thresholds = FALSE,
         nCharNodes = 0,
         mar = c(6, 6, 6, 6),
         color = list(lat = "lightblue", man = "white"),
         edge.color = "black",
         edge.width = 1.5,
         title = FALSE)
```

## 5. Report Fit Indices and Coefficients

```{r report-fit}
fitMeasures(fit, c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr")) %>%
  enframe(name = "Index", value = "Value") %>%
  mutate(Threshold = case_when(
    Index == "cfi" ~ ">= 0.95",
    Index == "tli" ~ ">= 0.95",
    Index == "rmsea" ~ "<= 0.06",
    Index == "srmr" ~ "<= 0.08",
    Index == "chisq" ~ "n.s. ideally",
    Index == "df" ~ "-",
    Index == "pvalue" ~ "> .05 ideally",
    TRUE ~ "")
  ) %>%
  kable(digits = 3, caption = "Table X. Model Fit Indices") %>%
  kable_styling(full_width = FALSE)
```

```{r report-coefficients}
parameterEstimates(fit, standardized = TRUE) %>%
  filter(op == "~") %>%
  select(lhs, op, rhs, est, pvalue, std.all) %>%
  rename(Standardized = std.all) %>%
  kable(digits = 3, caption = "Table Z. Standardized Path Coefficients") %>%
  kable_styling(full_width = FALSE)
```

## 6. Standardized Factor Loadings

```{r factor-loadings}
parameterEstimates(fit, standardized = TRUE) %>%
  filter(op == "=~") %>%
  select(lhs, rhs, std.all, se, pvalue) %>%
  rename(Latent = lhs,
         Indicator = rhs,
         `Std. Loading` = std.all,
         `SE` = se,
         `p-value` = pvalue) %>%
  kable(digits = 3, caption = "Table Y. Standardized Factor Loadings") %>%
  kable_styling(full_width = FALSE)
```

## 7. Composite Reliability & AVE (semTools)

```{r reliability-ave}
cr_vec <- compRelSEM(fit)
ave_vec <- AVE(fit)

df_cr <- data.frame(
  `Latent Variable` = names(cr_vec),
  `Composite Reliability` = round(as.numeric(cr_vec), 3),
  `Average Variance Extracted (AVE)` = round(as.numeric(ave_vec), 3)
)

kable(df_cr, caption = "Table AA. Composite Reliability and AVE") %>%
  kable_styling(full_width = FALSE)
```

## 8. Latent Correlation Matrix

```{r latent-correlation-matrix}
inspect(fit, "cor.lv") %>%
  round(3) %>%
  kable(caption = "Table BB. Latent Variable Correlation Matrix") %>%
  kable_styling(full_width = FALSE)
```

## 9. Discriminant Validity (Fornell–Larcker Criterion)

```{r discriminant-validity}
ave <- AVE(fit)
sqrt_ave <- sqrt(ave)
cor_lv <- inspect(fit, "cor.lv")
fornell <- cor_lv
for (i in 1:nrow(fornell)) fornell[i, i] <- sqrt_ave[i]
round(fornell, 3) %>%
  kable(caption = "Table CC. Fornell–Larcker Criterion (Diagonal = sqrt(AVE))") %>%
  kable_styling(full_width = FALSE)
```
